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Abstract. We present first results of the non-linear evolution of rotating relativistic 
stars obtained with an axisymmetric relativistic hydrodynamics code in a fixed 
spacetime. As initial data we use stationary axisymmetric and perturbed 
configurations. We find that, in order to prevent (numerical) angular momentum 
loss at the surface layers of the star a high-resolution grid (or a numerical scheme 
that retains high order at local extrema) is needed. For non-rotating stars, we 
compute frequencies of radial and non-radial small-amplitude oscillations, which are 
in excellent agreement with linear normal mode frequencies computed in the Cowling 
approximation. As a first application of our code, quasi-radial modes of rapidly 
rotating relativistic stars are computed. By generalizing our numerical code to 3-D, 
we plan to study the evolution and non-linear dynamics of toroidal oscillations (r- 
modes) of rapidly rotating neutron stars, which are a promising source of gravitational 
waves. 



1. Introduction 

The numerical evolution of neutron stars in full General Relativity has been the 
focus of many research groups in recent years [0, ||, ^, |^. So far, these studies have 
been limited to initially non-rotating stars. However, the numerical investigation of 
many interesting astrophysical applications, such as the rotational evolution of proto- 
neutron stars and merged neutron stars or the simulation of gravitational radiation 
from unstable pulsation modes, requires the ability of accurate long-term evolutions 
of rapidly rotating stars. We thus present here the first study of hydrodynamical 
evolutions of rotating neutron stars in the approximation of a static spacetime. This 
approximation allows us to evolve relativistic matter for a much longer time than 
present coupled spacetime plus hydrodynamical evolution codes. Since the pulsations 
of neutron stars are mainly a hydrodynamical process, the exclusion of the spacetime 
dynamics has only a limited effect and allows for qualitative conclusions to be drawn. 

The rotational evolution of neutron stars can be affected by several instabilities 
(see Q for a recent review). If hot protoneutron stars are rapidly rotating, they can 
undergo a dynamical bar- mode instability p]. When the neutron star has cooled to 
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about 10^°K after its formation, it can be subject to the Chandrasekhar-Friedman- 
Schutz instability j^, ^ and it becomes an important source of gravitational waves. 
It was recently found that the I = m r-mode has the shortest growth time of the 
instability |l^ and it can transform a rapidly rotating newly-born neutron star 
to a Crab-like slowly-rotating pulsar within about a year after its formation jll] , p^ . 
In this model, there are two important questions still to be answered What is 
the maximum amplitude that an unstable r-mode can reach (limited by nonlinear 
saturation) and is there any transfer of energy to other stable or unstable modes via 
non-linear couplings? Such questions cannot be answered by computations of normal 
modes of the linearized pulsation equations, but require non-linear effects to be taken 
into account. We therefore need to develop the capability of full non- linear numerical 
evolutions of rotating stars in General Relativity. 

Our present 2-D (axisymmetric) code uses high-resolution shock-capturing 
(HRSC) finite-difference schemes for the numerical integration of the general 
relativistic hydrodynamic equations |l4| (see fl^ for a recent review of applications 
of HRSC schemes in relativistic hydrodynamics). In a similar context to the one 
presented here let us note that such schemes have been succesfuUy used before in the 
study of the numerical evolution and gravitational collapse of non-rotating neutron 
stars in 1-D ||l6|. An alternative approach, based on pseudospectral methods, has 
been presented in [ p7| . 

Using our code in 1-D time-evolutions we can accurately identify specific normal 
modes of pulsation. In 2-D the code is suitable for the evolution of rotating stars, 
with the additional complication of having to pay special attention to an angular 
momentum-loss at the (non-spherical) surface of the star, as we will show below. 



2. Initial Configurations 

Our initial models are fully relativistic, stationary and axisymmetric configurations, 
rotating with uniform angular velocity Q. The metric in quasi-isotropic coordinates 
El is 

= -e^'^dt^ + B^e-^"r^ sin^ 9{d(j) - uodtf + e^°'{dr'^ + r^dO^), (1) 

where v, B, a and uj are metric functions (gravitational units are implied). In the 
non-rotating limit the above metric reduces to the metric of spherical relativistic stars 
in isotropic coordinates. 

We assume a perfect fluid, zero-temperature equation of state (EOS), for which the 
energy density is a function of pressure only. The following relativistic generalization 
of the Newtonian polytropic EOS is chosen: 

P^Kpa (2) 
e=po + Np, (3) 

where p is the pressure, e is the energy density, po is the rest-mass density, K is the 
polytropic constant and N is the polytropic exponent. 

The initial equilibrium models are computed using a numerical code by Stergioulas 
& Friedman Jl^ which follows the Komatsu, Eriguchi & Hatchisu |20j method (as 
modified in pl|) with some changes for improved accuracy (see [ p2[ for a comparison 
with other codes). The code is freely available and can be downloaded from the 



following URL address: http://www.gravity.phys.uwm.edu/Codc/rns 
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3. Relativistic Hydrodynamic Equations 

The equations of (ideal) relativistic hydrodynamics are obtained from the local 
conservation laws of density current, J** and stress-energy, T'^" 

V^J'' =0 (4) 

V^T^'^ = (5) 

with 

^P^u^^ (6) 

T^"^ = po/^^i^"^+TO''^ (7) 

for a general EOS p = p{p, e). This choice of the stress-energy tensor limits our study 
to perfect fluids. 

In the previous expressions is the covariant derivative, is the fluid 4- velocity 
and h is the specific enthalpy 

h = l + e + — (8) 
Po 

with £ being the specific internal energy, related to the energy density e by 

e^^-l. (9) 
Po 

With an appropriate choice of matter fields the equations of relativistic 
hydrodynamics constitute a (non-strictly) hyperbolic system and can be written 
in a flux conservative form, as was first shown in [ p3[ for the one-dimensional 
case. The knowledge of the characteristic fields of the system allows the numerical 
integration to be performed by means of advanced high-resolution shock-capturing 
(HRSC) schemes, using approximate Riemann solvers (Godunov-type methods). The 
multidimensional case was studied in [ p4| , within the framework of the 3-1-1 
formulation. Further extensions of this work to account for dynamical spacetimes, 
described by the full set of Einstein's non- vacuum equations, can be found in 
Fully covariant formulations of the hydrodynamic equations (i.e., not restricted to 
spacelike approaches) and also adapted to Godunov-type methods, are presented in 

11, H- 

In the present work we use the hydrodynamic equations as formulated in ||24| . 
Specializing for the metric given by Eq. (|l]), the 3-1-1 quantities read 

a =6" (10) 

= -ujB^e''^''r^sm'^9 (11) 

Irr = e^" (12) 

lee = r^e^" (13) 

70^ = B^-^^r^ sin' 9 (14) 

where a is the lapse function (the tilde is used to avoid confussion with the metric 
potential a) and /J^ is the azimuthal shift. 

The hydrodynamic equations are written as a first-order flux conservative system 
of the form 

du daf^ daf^ . 
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where u, P, f and s are, respectively, the state vector of evolved quantities, the radial 
and polar fluxes and the source terms. More precisely, they take the form 



U = {D, Sr, Sg, SiI),t) 

f = {Dv'-,Srv'- +p,Sev\S^v\{T+p)v'-) 
f = {Dv',Sy,Sgv''+p,S^v'',{T + p)v'). 

The source terms can be decomposed in the following way 
s = as —at at 



dr 



do 



with 7 = det 7ij and 



dg, 



(16) 
(17) 
(18) 



(19) 



(20) 



with j — r, 6, (j). The definitions of the evolved quantities in terms of the "primitive" 
variables (p, Vj , e) are 



D = poW 

Sj = pahW^Vj 

T = pohW^ -p- D 

where W is the relativistic Lorentz factor 
W = au^ 



1 



(21) 
(22) 
(23) 



(24) 



7^jV 



The 3-velocity components are obtained from the spatial 



components of the 4-velocity in the following way 

v' ^ h — . 

W & 



(25) 



Explicit expressions for the non- vanishing Christoffel symbols for metricO) , appearing 
in the source terms of the hydrodynamic equations, are presented in ||27|. 



4. Numerical Methods 

As stated before, our numerical integration of system (|l5|) is based on Godunov-type 
methods (also known as HRSC schemes). In a HRSC scheme, the knowledge of the 
characteristic fields (eigenvalues) of the equations, together with the corresponding 
eigenvectors, allows for accurate integrations, by means of either exact or approximate 
Riemann solvers, along the fluid characteristics. These solvers, which constitute 
the kernel of our numerical algorithm, compute, at every interface of the numerical 
grid, the solution of local Riemann problems (i.e., the simplest initial value problem 
with discontinuous initial data). Hence, HRSC schemes automatically guarantee 
that physical discontinuities appearing in the solution, e.g., shock waves, are treated 
consistently (the shock- capturing property). HRSC schemes are also known for giving 
stable and sharp discrete shock proflles. They have also a high order of accuracy, 
typically second order or more, in smooth regions of the solution. 
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We perform the time update of system ( [15D according to the foUowing conservative 
algorithm: 

= K,j - ^(fi+l/2j ~%-l/2.j) 

- ^(M^,3 + l/2 - g»j-l/2) + AtS.j . (26) 

Index n represents the time level and the time (space) discretization interval is 
indicated by At ( Ar, A6') . The "hat" in the fluxes is used to denote the so-called 
numerical fluxes which, in a HRSC scheme, are computed according to some generic 
flux-formula, of the following functional form (suppressing index j): 

\±\ - \ (uf±i) + f(uf±i) - E I V I Ai3„F„^ . (27) 

Notice that the numerical flux is computed at cell interfaces (i±I/2). Indices L and R 
indicate the left and right sides of a given interface. Quantities A, Aw and r denote the 
eigenvalues, the jump of the characteristic variables and the eigenvectors, respectively, 
computed at the cell interfaces according to some suitable average of the state vector 
variables. Generic expressions for the characteristic speeds and eigenfields can be 
found in [Q. Our code has the ability of using different approximate Riemann solvers: 
the Roe solver psf , widely employed in fluid dynamics simulations, with arithmetically 
averaged states and the Marquina solver (29), which has been extended to Relativity 
in |3C[| . The computations presented here were obtained using Marquina's scheme. 

A technical remark: the equilibrium star is supplemented by a low-density 
uniform atmosphere, which is necessary for computing non-singular solutions of the 
hydrodynamic equations everywhere in the computational domain. After each time- 
step we reset the atmosphere's density and pressure to their initial values, avoiding 
unwanted accretion of matter onto the star. The influence of the atmosphere is thus 
restricted to the surface grid-cells. 



5. Pulsations of Non-rotating Stars 

Since our code uses spherical polar coordinates, it can also be employed to study the 
evolution of non-rotating stars in 1-D. In the evolution of initially static non-rotating 
stars, we observe the following properties (note that our numerical grid is Eulerian): 

(i) Small-amplitude radial pulsations are triggered by the truncation errors of the 
finite-differencing scheme. 

(ii) The radial pulsations arc dominated by a set of discrete frequencies, which 
correspond to the normal modes of pulsation of the star. 

(iii) The numerical viscosity of the finite-difference scheme damps the pulsations and 
the damping is stronger for the higher frequency modes. 

(iv) The presence of a constant density atmosphere affects the finite differencing at 
the surface grid-cells, which increases the numerical damping of pulsations and 
also causes a continuous but very small drifting of the density distribution. 
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The initial amplitude of the radial pulsations and the small drift in density converge 
to zero at a second order rate with increasing resolution. The value of the density in 
the atmosphere region has a large effect on the damping of the pulsations. If it is too 
large, the damping is strong. To minimize this effect, we typically set the density of 
the atmosphere equal to 10"^ times the density of the last grid point inside the star. 

5.1. 1-D evolutions 

We study the numerical evolution of a nonrotating iV = 1.5 relativistic polytrope 
with M/R = 0.056. The star is immediately set into radial pulsation, triggered by 
the finite difference truncation errors. The time-evolution of the radial velocity Vr, 
computed at 25% of the star radius using a grid of 400 zones, is shown in Fig. |^. The 
vertical axis is dimensionless (c = G = A/q = 1). The radial velocity is initially a 
very complex function of time. As we will show next, the pulsation consists mainly 
of a superposition of normal modes of oscillation of the fluid. The high frequency 
normal modes are damped quickly and after 20ms the star pulsates mostly in lowest 
frequency modes. Because these oscillations are caused only by the truncation errors, 
the magnitude of the radial velocity is extremely small. As shown in Fig. |^ it is 
only a few times larger than the non-zero residual velocity around which the star is 
oscillating. This residual velocity converges to zero as second order with increased 
resolution. 
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Figure 1. Evolution of the radial velocity of an initially static, non-rotating star, 
caused by the truncation errors of the finite-difference scheme. The radial pulsations 
are mainly a superposition of normal modes of the star. 

The small-amplitude radial pulsations in the non-linear, fixed spacetime evolutions 
correspond to linear normal modes of pulsation in the relativistic Cowling 
approximation, in which perturbations of the spacetime are ignored. A Fourier 
transform of the density or radial velocity time-evolution can be used to identify the 
normal mode frequencies. Fig. shows the Fourier transform of the radial velocity 
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evolution shown in Fig. The normal mode frequencies stand out as sharp peaks 
on a continuous background. The width of the peaks increases with frequency. The 
frequencies of radial pulsations identified from Fig. are shown in Table 



0.013 



0.010 - 



0.007 



0.005 



0.003 - 



0.000 



HI: H2 H3 



H4 



H5i 



H6 H7 H8 H9 



0.0 2.0 4.0 6.0 8.0 10.0 
f (kHz) 



Figure 2. Fourier transform of the evolution of the radial velocity in Fig. gj. The 
frequencies are in excellent agreement with linear normal mode frequencies computed 
with an eigenvalue code. The units of the vertical axis are arbitrary. 

To compare the obtained frequencies to linear normal mode frequencies, we use a 
different code that solves the linearized relativistic pulsation equations for the stellar 
fluid, in the Cowling approximation |3|, |3^ , as an eigenvalue problem. In Table 
we present the results of this comparison. The typical agreement between frequencies 
computed by the two methods is better than 0.5% for the fundamental F-mode and the 
lowest frequency harmonics Hi — H-^ and better than 0.8% for the the higher harmonics 

— Hg . This is a strong test for the accuracy of the evolution code and our results 
can be used as a testbed computation for other relativistic multi-dimensional evolution 
codes. 

5.2. 2-D evolutions 

In a similar way, small-amplitude non-radial pulsations can be studied with the 
present evolution code and the obtained frequencies can be compared to perturbation 
results. We find that the truncation errors of the finite difference scheme do 
not excite non-radial pulsations to a sufficiently large amplitude compared to the 
amplitude of radial pulsations, so that one cannot identify them accurately in a 
Fourier transform. Instead, one has to perturb the initial configuration, using an 
appropriate eigenfunction for each nonradial angular index Such a perturbation can 
be constructed using the eigenfunctions of linear pulsation modes, computed with the 
perturbation code in the Cowling approximation. The frequencies of the non-radial 
modes are then found from a Fourier transform of the time-evolution of the velocity 
component ve. 
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Radial Pulsation Frequencies 



Mode non-linear code (kHz) Cowling (kHz) difference 



F 


1.703 


1.697 


0.3% 


Hi 


2.820 


2.807 


0.5% 


H2 


3.862 


3.868 


0.02% 


H3 


4.900 


4.910 


0.2% 


Hi 


5.917 


5.944 


0.4% 


H5 


6.930 


6.973 


0.6% 


He 


7.947 


8.001 


0.7% 


Hr 


8.960 


9.029 


0.8% 


Hs 


9.973 


10.057 


0.8% 


Hg 


11.030 


11.086 


0.5% 



Table 1. Comparison of small-amplitude radial pulsation frequencies obtained with 
the present non-linear evolution code to linear perturbation mode frequencies in the 
relativistic Cowling approximation. The equilibrium model is a nonrotating A'^ = 1.5 
relativistic polytrope with M/R = 0.056. 



Quadrupole Pulsation Frequencies 



Mode non-linear code (kHz) Cowling (kHz) difference 



/ 


1.28 


1.286 


0.5% 


Pi 


2.68 


2.681 


0.04% 


P2 


3.65 


3.699 


1.3% 


P3 


4.66 


4.719 


1.3% 


Pi 


5.66 


5.742 


1.4% 


P5 


6.83 


6.764 


1.0% 


P6 


7.80 


7.788 


0.2% 



Table 2. Comparison of small-amplitude quadrupole (Z = 2) pulsation frequencies 
obtained with the present non-linear evolution code to linear perturbation mode 
frequencies in the relativistic Cowling approximation. The equilibrium model is a 
nonrotating A' = 1.5 relativistic polytrope with M/R = 0.056. 



Table || shows a similar comparison as in Table |l| for the quadrupole {I = 2) 
pulsations of the same = 1.5 relativistic polytrope. Since the non-radial modes 
have to be computed on a 2-D grid, we cannot use resolutions as high as in the 1-D 
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computations. For a small grid-size of 80 x 80 zones and a total evolution time of 6ms, 
the agreement between frequencies computed by the two methods is better than 1.4% 
for the fundamental /-mode and the p-modes pi — pe- For this grid-size, frequencies 
higher than the p^ mode could not be computed accurately, because the grid is to 
coarse to resolve their eigenfunctions (higher harmonic eigenfunctions have a larger 
number of nodes in the radial direction). 

6. Rotating Stars 

We now turn to the evolution of initially stationary, uniformly rotating neutron stars. 
In these evolutions, we observe the same qualitative properties as for non-rotating 
stars (section ||) and an additional important property: the angular momentum of the 
star is not conserved at the surface layer. This is due to the fact that the velocity 
component of the fluid has a maximum at the surface, while the numerical scheme 
(although second-order accurate in smooth regions of the solution) is only first-order 
accurate at local extrema. Moreover, the code evolves the relativistic momenta. Si, 
and the velocity components (as well as the rest of "primitive" variables) must be 
recovered through a root finding procedure which involves dividing by the density. At 
the surface of the star (where the density is very small) this contributes to obtaining 
less than second-order accuracy. 



1.5 




r (km) 

Figure 3. Time-evolution of the velocity component of a rotating star (see text 
for details). The accuracy is second order in the interior but only first order at the 
surface. This results in an angular momentum loss of the surface layers. 

A representative example of the evolution of a rotating star is presented in Fig. 
^, which shows the evolution, at different times, of the velocity component w^. The 
star is again a = 1.5 polytrope with the same central density as the non-rotating 
star presented in section |^ and rotating at 74% of the mass-shedding limit at same 
central density. The evolution was for one rotation period on a 96 x 60 grid. The 
vertical axis is dimensionless {c — G — Mq — 1). The figure shows that the </>- velocity 
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Quasi-radial Pulsation Frequencies 



n/^K F (kHz) HI (kHz) H2 (kHz) 



0.0 


1.71 


2.82 


3.90 


0.32 


1.70 


2.77 


3.81 


0.44 


1.67 


2.71 


3.68 


0.62 


1.64 


2.52 


3.46 


0.74 


1.53 


2.38 


3.33 


0.84 


1.36 


2.25 


3.17 



Table 3. Frequencies of quasi-radial modes for different values of the ratio of angular 
velocity of the star Q, to the angular velocity at the mass-shedding limit VIk, for 
a sequence of rotating relativistic stars of same central density. The non-rotating 
member of the sequence is the same as in Table ll 



in the interior of the star remains close to its initial value, while it decreases as a 
function of time in the outer layers. We find that this is a generic property of the 
present numerical scheme for any rotation rate and for any grid-size. By comparing 
evolutions with different grid sizes, we verified that the loss of angular momentum at 
the surface improves as first-order with resolution, while the evolution of the velocity 
in the interior is second-order accurate. However, as the evolution proceeds in time, 
the first-order surface effect gradually affects the interior of the star. 

7. Quasi-radial Modes of Rotating Stars 

As a first application of our code, we compute quasi-radial modes (i.e. modes that 
in the non-rotating limit reduce to radial modes) of rapidly rotating relativistic stars 
in the Cowling approximation. Previously, these modes have been computed for fully 
relativistic stars only in the slow-rotation limit (but without the assumption of a fixed 
spacetime) by Hartle & Friedman jsi) (see also Q). We compute the three lowest- 
frequency quasi-radial modes for a sequence of rotating stars of same central density. 
The non-rotating member of the sequence is the non-rotating star of section |^. Table 
^ and Fig. ^ show our results for a low resolution grid of 100 x 80 zones (note that 
our computational grid assumes equatorial plane symmetry). For this resolution we 
estimate the accuracy of the frequencies to be of the order of 1 — 2%. 

For the sequence of stars considered here, the frequencies of the quasi-radial 
modes decrease with increasing rotation rate. This agrees with previous slow-rotation 
computations which predict a decrease as Vfi , where is the angular velocity of the 
star. For fast rotation, the change in the frequencies of quasi-radial modes is affected 
by higher order terms in 17, because of the large deformation of the equilibrium star. 
Also, for rapidly rotating stars the quasi-radial mode frequencies are more "closely 
packed" than in non-rotating stars. 
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Figure 4. Frequencies of tlie lowest three quasi-radial modes vs. the ratio of angular 
velocity of the star H to the angular velocity at the mass-shedding limit , for the 
sequence of rotating relativistic stars in Tabic M. 



8. Discussion 



Our axisymmetric relativistic hydrodynamical code is capable to evolve rapidly 
rotating stars in a fixed spacetime. We find that, for non-rotating stars, small 
amplitude oscillations have frequencies that agree with linear normal mode frequencies 
in the Cowling approximation and we compute the quasi-radial modes of rapidly 
rotating stars. 

Modern HRSC numerical schemes (as the ones used in our code), satisfying 
the "total variation diminishing" (TVD) property ]36| , are second-order accurate in 
smooth regions of the flow, but only first-order accurate at local extrema. In our 
rotating stars runs we find that this results in a loss of angular momentum of the 
surface layers of the star, which gradually also affects the interior of the star. This 
angular momentum loss only vanishes as first-order with incresing resolution and we 
thus conclude that for accurate long-term evolutions of rotating neutron stars it is 
essential to use rather fine grids. Furthermore, to reduce the computational cost, 
one could use surface-adapted coordinates or fixed-mesh refinement. It would also 
be interesting to see whether the loss of angular momentum per rotation period 
will be significantly smaller in a frame co-rotating with the star. An alternative 
solution to this problem, which we plan to investigate, could be the use of "essentially 
non-oscillatory" (ENO) schemes, which maintain high-order of accuracy even at local 
extrema pTj . 

All previous considerations are important for the study of the non-linear dynamics 
of unstable toroidal oscillations (r-modes) in 3-D, which have a long growth time and 
thus require highly accurate long-term evolutions. 
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